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Abstract 



^ ■ We determine the parameter Csw required for 0(a)-improvement of the three 

flavor Wilson fermion action together with the tree-level Symanzik improved gauge 

^ . action. The standard improvement condition is employed for a range of couplings. 

Additionally, we perform a check of the volume independence of Cgw and provide 

Q^i a preliminary estimate of the lattice spacing at our largest values of g^. 
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1 Introduction 

The continuum limit is an essential part of lattice QCD calculations. In order to 
control this extrapolation in the lattice spacing a, it is advisable to work with a 
discretization whose leading cut-off effects are O(a^). The theory and reduction of 
scaling violations in on-shell quantities is well-established p!P][5] and for Wilson 
fermions [4J 0(a) improvement can be achieved by adding a single dimension- 
five operator to the action [5]. This requires the non-perturbative tuning of its 
coefficient Csw and can be performed using a standard procedure Pf7] . Simulations 
with the resulting action in two-flavor flavor QCD have exhibited the expected 
moderate scaling violations [S]. 

Here we present a determination of the coefficient Cgw using the standard 
procedure for N{ = 3 flavor QCD with the tree-level improved Liischer-Weisz 
gauge action p]. We choose this gauge action because it has been demonstrated 
to possess superior scaling properties in pure gauge theory [TO] . 

This paper is organized as follows. After deflning our lattice regularization in 
Sec. m we summarize the standard improvement programme in Sec. [21 In Sec. [H 
details of the numerical computations are given together with the resultant inter- 
polating formula for the improvement coefficient Csw(fi'o)- We also give results for 
a renormalized quantity which suggests that the range of bare couplings covered 
by our simulations extends at least to a lattice spacing of a ~ 0.09 fm. 

2 Lattice setup 

The 0(a)-improved Wilson Dirac operator [1[[5] is given by 

-.3 3 . 

^w = 2 Y.^1,^K + ^A.) - ^^^m) + Csw Yl i^^-^M- + "^0 (2.1) 

with V^ and V* the covariant forward and backward derivatives, respectively, and 
F^,^ the standard discretization of the field strength tensor [llj. The bare mass 
niQ will be replaced below by the hopping parameter k, with tuq = {k,^^ — 8)/2. 

The gauge action Sq contains sums over all oriented 1x1 plaquettes as well 
as all 1 X 2 rectangles, which are denoted by the sets Sq and Si, respectively, 

5g = /3 5^ Cfc ^ WkiC) tr{l - U{C)}, (2.2) 

A:=0,1 CeSk 

where /3 = G/Qq, Cq = 5/3 and ci = —1/12. The weight factor Wk{C) is set to unity 
for all loops away from the boundaries. Schrodinger functional [T21IT3] boundary 
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conditions are imposed on the gauge fields such that 

, ^, all links in C are on a boundary 
MC) = \'' ^^ . (2.3) 

I 1, otherwise 

and 



1 

2' 



all links in C are on a boundary 
u7i(C) = ■^ |, C has exactly two links on a boundary (2.4) 

1, otherwise. 

This corresponds to 'Choice B' of Ref. [H] and guarantees boundary 0(a) im- 
provement at tree-level of perturbation theory. We also set the fermionic boundary 
counter term according to the 1-loop formula 



cf = 1- 0.0122 Cf gl with Cp = 4/3 . (2.5) 

Finally, the values of the spatial links at the boundaries are fixed to 

U{x, k) 1^.0=0 = expjaCfc} , (7^ = _- diag(-7r, 0, vr) , (2.6) 

t/(x, A;)|^(,=r = exp{aC^} , C^ = — diag(-57r, 27r, 37r) , (2.7) 

while the fermion fields satisfy periodic boundary conditions in the spatial direc- 
tions and the standard Schrodinger functional boundary conditions in time. 



3 Improvement condition 

The standard 0(a) improvement programme relies on the PCAC relation, which 
involves the improved axial-vector current (^i)" and the pseudoscalar density P" 
given by 

[A,)l = Al + ac^\{d; + d,)P\ (3.1) 



— i){x), P\X) = ^(X)75y 



A;:(x) = t/j{x)j,j,-ij{x), P\x) = ^(x)75-^(x) . (3.2) 



In the unimproved theory the unrenormalized PCAC relation 

hd, + d;){{Ai)l{x)0) = 2m{P'^{x)0) (3.3) 



is violated by terms of 0(a). By using three different choices oi x, O and O' one 
can define Cgw and ca, requiring that m is the same in all three cases. In this 
situation Eq. 13. f I holds up to O(a^). 

This method has been applied to the Nf = 0,2,3,4 cases for a variety of 
different actions Pl7|15fl6fl7fl8] . For technical reasons, we use lattices with T = 
2L — a. However, this additional 0(a) effect is irrelevant, as the determination 
of Csw using our improvement condition is ambiguous at 0(a). We employ the 
PCAC relation with boundary operators O and O' on time slices Xq = and 
xq = T, respectively, 

0^ = a'Yl C(y)75y C(z) , O" = a'Yl C'(y)75y C'(z) • (3.4) 



The correlation functions which enter the PCAC relation Eq. 13.11 are then 

/a(xo) = -^(AS(x)0'^) , /p(xo) = -^(P'^(x)0«) , (3.5) 

/;(r - xo) = +i(Ag(x)0"^) , /^(T - xo) = -i(P«(x)0'") , (3.6) 

As has been suggested in Ref. [6], effective masses M(xo) defined by 

M(xo, yo) = r(xo) - sixo f}/^']'"^^^'] . (3.7) 

s'iyo) - s{yo) 

with 

rixo) = -{Oq + do)fAixo)/fp{xo) and s(xo) = -a^o^o /p(xo)//p(xo) (3.8) 

correspond to a particular choice of ca in the improved currents. Since M{xo,yo) 
renormalizes multiplicatively, it is a useful quantity to define the quark mass at 
fixed P and Cgyf. Specifically, we take as our definition of the quark mass 

M = ^ {M{L, L/2) + M{L - a, L/2)) . (3.9) 

Tuning M k. {] defines the values of k at which we impose the improvement 
condition. 

As discussed above, this improvement condition requires the difference be- 
tween two masses to vanish. In addition to M(xo,i/o), a second mass M'(xo,i/o) 
is considered where r and s in Eq. 13.71 are replaced by their primed counterparts. 
These masses are evaluated at Xq = 3T/4 and yo = T/4. However, because we 
have T = 2L — a, we round the two arguments towards the center of the lattice. 



The parameter Cgw is chosen such that the difference between these two masses 
AM is equal to its tree-level value AM*^°) 

AM = M(3T/4, T/4) - M'(3T/4, T/4) = AM(°) . (3.10) 

The numerical value of AM*^°) depends on the lattice geometry and can be 
computed using the solution to the Dirac equation as given in Sec. 6.2 of Ref. [T^ . 
For the 15 x 8^ lattices aAM(°) = 0.000393. Additionally, this offset may be 
obtained from measurements on free gauge fields. A stringent test of our entire 
workflow is the reproduction of the AM'^''^ obtained from analytic calculations 
using simulations at large values of /3. 

In principle it would be preferable to keep the physical size of the system L 
constant as the continuum limit is approached. This, however, turns out to be very 
costly in practice, because autocorrelations associated with the topological charge 
sectors quickly become very large as the lattice spacing is lowered. We therefore 
opt to impose the improvement condition at fixed L/a, where the contribution 
of non-zero topological charge sectors decreases rapidly in the continuum limit. 
In Sec. 15.31 we will show that the effect of the finite volume does not seem to be 
relevant at the current level of accuracy. 

4 Simulations 

For the simulations we use the openQCD code, which is publicly availible onlinci 
and implements the lattice setup of Sec. |2]for several simulation algorithms. These 
algorithms are the subject of Ref. [22] so we restrict ourselves to a brief summary 
here. 

We employ the HMC algorithm [2lJ with a twisted-mass Hasenbusch fre- 
quency splitting [22|23] for a doublet of two of the three degenerate quarks. For 
most ensembles with /3 < 3.5 twisted mass reweighting [21] is used, i.e. we simulate 
with a small twisted mass /i = 0.001 and then include a stochastically estimated 
reweighting factor to correct for this in the measurement. This significantly in- 
creases the stability of the simulation. The third quark is simulated using the 
RHMC algorithm [25|. In all cases, we use a nine pole rational approximation in 
the interval [0.02, 7.2] and correct for the rational approximation with an addi- 
tional stochastically estimated reweighting factor. All fermion determinants are 
factorized using even-odd preconditioning. 

For the entire range of couplings, we generated 15 x 8^ lattices at a variety of 
Csw values listed in Tabled] In addition, 23 x 12^ ensembles at /3 = 3.8 have been 
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generated as finite volume checks. To obtain a preliminary estimate of tlie lattice 
spacing at /? = 3.3 and 3.4, we performed some L = T runs (using a modified 
version of the openQCD code) with L/a = 8 and zero boundary fields. These 
runs used both our discretization and the one of Ref. [26] and are summarized in 
Table [31 

For the L/a = 8 lattices at the four smallest values of /3 we use a three-level 
hierarchical integration scheme in which the outermost level employs the second 
order Omelyan-Mryglod-Folk (OMF) integrator, while the two inner levels use the 
fourth order OMF integrator [27] . The force from the pole closest to the origin in 
the rational approximation is integrated on the coarsest timescale, the remaining 
fermion forces on the intermediate timescale, and the gauge force on the finest 
timescale. Four or five outermost iterations together with a single iteration of 
the remaining two integrators typically achieve ~ 90% acceptance for molecular 
dynamics trajectories of length r = 2. 

For the L/a = 12 runs at /3 = 3.8, the remaining L/a = 8 runs and the T = L 
runs, a two level scheme with both levels using fourth order OMF was found to 
be effective, with between 5 and 8 steps in the outermost integrator and a single 
inner iteration. 

The integration schemes discussed above were very stable in all cases, result- 
ing in small Hamiltonian violations. The most difficult simulations were those 
with L/a = S,T = 2L — a at (3 = 3.3 at the smallest value of Cgw, but even those 
had only about 0.6% trajectories with AH > 10. At Cgw = 2.1 this falls already 
to 0.2% and we had 0.08% of such events at Cgw = 2.4. At the smallest Cgw = 1-7 
for P = 3.4 these trajectories occurred with 0.1%, a percentage rapidly falling for 
larger values of Csw and (3. 

5 Results 

We finally come to the results of our simulations. Our analysis strategy is discussed 
in Sec. 15.11 and results for Csw(fl'o) are collected in Sec. 15.21 with a finite volume 
check in Sec. 15.31 A preliminary scale determination is given in Sec. 15.41 

A summary of the ensembles generated for the determination of Csw appears 
in Tab. [H where we also give the total statistics accumulated over several replica. 
These ensembles consist of L/a = 8 lattices used for our final result as well as those 
used for the finite volume check. The ensembles generated for the preliminary scale 
setting are discussed in Sec. 15.41 and collected in Tab. [31 
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-0.0092(16) 


0.0028(8) 
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8 '-'^ 
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-0.0010(12) 
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0.1270442 


0.0070(3) 


-0.0050(2) 
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0.1374470 


0.0091(7) 


0.0004(4) 


48860 


3.5 


8 2.20 


0.1319060 


0.0010(3) 


-0.0020(2) 


62872 




2.55 


0.1267940 


0.0109(4) 


-0.0056(3) 


37800 




1.50 


0.1420500 


-0.0009(14) 


0.0030(7) 


19200 


■^ fini 


8 '-^^ 


0.1387200 


-0.0057(9) 


0.0006(5) 


22200 




^ 1.90 


0.1353560 


-0.0010(6) 


-0.0002(4) 


20200 




2.10 


0.1319920 


0.0130(4) 


-0.0027(3) 


28600 




1.20 


0.1434300 


0.0134(10) 


0.0039(5) 


47400 


s s 


8 '-'^ 


0.1405500 


0.0041(7) 


0.0023(5) 


21942 




^ 1.60 


0.1376520 


0.0006(5) 


0.0011(3) 


25762 
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0.1319400 


0.0043(4) 


-0.0029(2) 


43986 




1.20 


0.1434300 


0.0095(5) 


0.0012(3) 


14967 


3.8 


12 1.60 


0.1376520 


0.0005(3) 


0.0003(2) 


7000 




2.00 


0.1319220 


0.00952(15) 


-0.00095(16) 


12500 




1.00 


0.1410350 


0.0071(13) 


0.0045(5) 
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0.1374570 


0.0053(5) 
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0.00517(16) 


0.00053(14) 
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1.40 


0.1310350 


0.00394(17) 


-0.00174(13) 
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1.60 


0.1296320 


0.00118(11) 


-0.00397(11) 


7174 



Table 1: Simulation parameters for the runs used in the Csw determination as 
well as the resultant values for aM and aAM calculated on those ensembles. The 
integrated molecular dynamics time of all replica for each ensemble is also given. 
For all of these runs the boundary fields are as specified in Eq. 



P 3.3 3.4 3.5 3.601 3.8 4.3 6.0 

csw 2.13(5) 1.96(4) 1.90(3) 1.78(3) 1.61(2) 1.43(2) 1.213(9) 

Table 2: Values of the optimal Csw for which aAM = aAM^^\ Each value is from 
a linear interpolation of L/a = 8 data at a fixed value of /3. 

5.1 Analysis 

For each value of the coupling, we have to find the value of Cg^ and the quark 
mass for which the improvement condition is satisfied. It has already been shown 
in the past that the condition M = does not have to be achieved to a very high 
accuracy [7|T7] and we therefore require \aM\ < 0.015 for each individual point. 
We then measure aAM for several choices of Cgw and interpolate linearly to find 
the point with AM = AM^^\ 

The measurements of the required fermionic correlation functions (Eq. 13.51) 
are separated by one trajectory of length r = 2. Using the methods and software 
of Ref. [28j, we determine the integrated autocorrelation times of the observables 
aM and aAM and find them to be at most around 8 units of molecular dynamics 
time such that all ensembles provide at least 500 independent measurements and 
thus a reliable determination of the errors. This is confirmed by the normal 
distribution of mean values on single replica. 

Also, by studying the observables defined through the Wilson flow [29] — 
the action density constructed from smoothed links and the topological charge — 
we ensure that field space is sufficiently sampled in all simulations. In particular 
the topological charge is known to cause problems in the continuum limit [30]. 
However, since we keep L/a fixed instead of L, the contribution of sectors of 
non-zero topological charge diminishes rapidly as a — )■ 0. 

The integrated autocorrelation times of the smoothed action and the topo- 
logical charge at various smoothing ranges do not exceed 100 units of molecular 
dynamic time in any of our simulations. We therefore conclude that also for these 
observables configuration space is sampled sufficiently. 

5.2 Csw for L/a = 8 

The results for the improvement condition as a function of Cgw as well as the 
resultant linear interpolations are shown in Fig. [T]for six values of /3. The level 
of statistics is such that Cgw is determined with better than 3% precision in all 
cases. The optimal Cgw values at fixed (3 are collected in Tab. |2] and shown in 
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Figure 1: The improvement condition AM for several values of Cgw at each value 
of p. Linear fits to the improvement condition as a function of Cgw at fixed /3 
are also shown. The resultant values of Csw which minimize the improvement 
condition are given in Table [2l 



Fig. |2] together with the best-fit interpolating curve 

1 - 0.1921 g^ - 0.1378 c/^ + 0.0717^^ 



--swVi/o 



(g'o 



1 - 0.3881 gi 



(5.1) 



We see that the fit describes the data well and that the data approaches the known 
1-loop result at large (3, which we use to constrain the interpolating curve. 

5.3 Finite volume check 

Although Csw has typically been determined at fixed L/a = 8, it is interesting 
to assess the impact of this choice on the final result. As stated above, it would 
actually be preferable to keep L fixed when imposing the improvement condition. 
Smaller volumes are advantageous as they are computationally cheaper and have 
a larger slope in AM vs. Cgw, resulting in a clearer signal. Also, the problem 
of autocorrelations associated with freezing topological modes is much reduced. 
However, small volumes give rise to potentially large 0(a) effects in Cgw 

We detail here a single check at /3 = 3.8 between L/a = 8 and L/a = 12. The 
result is shown in Fig. [31 Within the statistical accuracies, the two volumes give 
the same value of Csw With Csw = 1-64(2) for L/a = 8 and 1.61(5) for L/a = 12, 
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Figure 2: Calculated values for Cgw together with the interpolating function repre- 
sented by the solid line. The dashed line is given by one-loop perturbation theory. 
To quantify finite volume effects, a value from simulations at 23 x 12'^ is given, 
with the value for Qq slightly shifted to the left for clarity. 



one can conclude that at least up to that lattice spacing, the smaller lattices are 
a good choice for the determination of Cqw The point at L/a = 12 also does not 
deviate by more than one standard deviation from the interpolating curve plotted 
in Fig. m Therefore, given our statistical accuracy, the systematic effects from 
fixing L/a = 8 do not appear to be significant. 



5.4 Preliminary scale determination 

To estimate the lattice spacing we compute a quantity with our discretization as 
well as the discretization used by PACS-CS, where the scale is known from large 
volume simulations using the mass of the Q baryon [31j. For comparison we use 
the couphng defined in Ref . |22] , except with periodic spatial boundary conditions 
for the fermion fields, i.e. with 6 = 0. This coupling is a renormalized quantity 
and at M = depends only on L, up to scaling violations. This lattice size L at a 
given value of the coupling serves as the dimensionful quantity to set the physical 
scale. 

Our results for this coupling are tabulated in Table O In both discretizations 
we expect cutoff effects of order 0{agQ) as boundary improvement for the gauge 
fields is implemented at tree level only. Furthermore, boundary improvement for 
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Figure 3: A comparison of the improvement condition at L/a = 8 and L/a = 12 
for (3 = 3.8. The values of Cgw are shghtly displaced for better clarity. 



the fermion fields is implemented only at 1-loop, resulting in additional Olag^) 
effects. We see that the f3 = 3.3 result for the coupling in our discretization lies 
above the value at a = 0.09fm, while the f3 = 3.4 is below. This suggests that the 
P corresponding to a = 0.09fm in our discretization is in the range [3.3, 3.4]. 



6 Conclusions 

In this paper we have determined Cgw for Nf = 3 lattice QCD with the tree- 
level Symanzik-improved gauge action. The result of our determination is the 
interpolation formula 



c 



{9l 



1 - 0.1921 gl - 0.1378^^ + 0.0717 1/^ 



swvyo 



(6.1) 



1-0.3881^2 

which may be taken as a definition of the lattice action. While our determination 
was performed at fixed L/a = 8, we performed a finite volume check at /3 = 3.8 
and found no significant change in Csw 

In addition to a determination of Csw, we have calculated a renormalized in- 
dependent coupling at /9 = 3.3 and 3.4 and compared with an alternative Nf = 3 
discretization. This determination suggests that the bare coupling correspond- 
ing to a ^ 0.09fm in our discretization is located in the interval (3 G [3.3,3.4], 
indicating that our Csw determination spans the range of desired lattice spacings. 
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Csw 


Cp 
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a(fni) 


9\L) 


MDU 
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1.9 


-0.331 


1.715 


0.972168 
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2.127114 
1.986246 


0.970424 
0.971294 


0.137017 
0.137553 


— 


7.381(72) 
6.225(21) 


21904 
38408 



Table 3: Nf = 3 results for the Wilson flow coupling g'^{L) from Ref. [52] using 
the Iwasaki gauge action and our discretization. We take the parameters and the 
scale determination from Ref. [31], while k was taken from Ref. [20] and found 
to give aM ^ 0. For these simulations only, we set T = L with boundary fields 
Ck = C', = 0. 
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